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ABSTRACT 


Ekman's linear equations for time dependent flow (neglecting 
wind stress) are solved using a time dependent Green's function and 
the method suggested by Welander (1957). The solution represents 
the vertical velocity profile in terms of the local time history of 
the changes in sea surface elevation determined by the divergence 
of the flow in the vertically integrated continuity equation. 

A fully implicit finite-difference scheme is developed to re- 
present a time dependent seiche oscillating across a shallow infinite 
channel. The transients associated with the formation of the bottom 
Spiral are clearly represented by the model and the influence of 
friction and Coriolis are individually and collectively introduced. 
The model allows independent calculation of velocity, volume trans- 
port, sea surface elevation, bottom stress, and the total energy 
balance of the system. The numerical scheme provides a method for 
adequately describing and investigating certain classes of time 
dependent motion, and its development suggests a mechanism for 


improving numerical prediction of storm surge. 
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I. INTRODUCTION 


"The purpose of computing is insight not numbers" 
(Hamming, 1962) 

A. BACKGROUND 

Oceanographers use a number of techniques for the study of 
oceanic flow including analysis of hydrographic data, use of 
hydraulic laboratory models, and mathematical treatment of the 
basic equations of motion. Initially oceanic circulation was 
treated as a steady state problem, and various investigators have 
been successful in representing and explaining many of the gross 
features of oceanic circulation using each of these techniques. 

Mathematical treatment of oceanic circulation began with 
Ekman's analysis (1905) of wind-driven circulation and his results 
have been extended by a number of authors. (Sverdrup, 1947, Munk, 
1950, Bryan, 1959) For the most part the study of large scale 
ocean dynamics has been done with steady state formulations. 

Mathematical treatment of time dependent motion is far less 
complete in the literature than that of steady state. It is also 
far more difficult to compare results due to the scarcity of 
reliable data. The complex nature of the solutions to time de- 
pendent boundary value problems further increases this difficulty. 
These solutions take the form of varying differential, integral, 
and infinite series equations whose numerical evaluation was not 
feasible until the advent of the high speed computer and the de- 


velopment of appropriate numerical techniques. 





In the last decade these methods have been applied, through 
numerical models, to the problems of explaining and predicting the 
various details of fluid motion in both the atmosphere and the 
oceans. Although much success has been attained with these methods, 
many critical questions still need explanation. Many of these 
questions relate to the structure of the various frictional boundary 
layers, the details of their development, and the resulting effects 
on the interior flow. 

Many complicated engineering problems occur in near shore and 
other shallow water regions. The time dependent nature of these 
areas is well documented and the resulting flow is greatly modified 
due to interaction with the ocean bottom. Therefore an understanding 
of the nature of time dependent flow at all levels is fundamental to 
dealing with problems in these areas. The question of how readjust- 
ment takes place and how the many transients can be appropriately 
represented in the bottom frictional boundary layer has not been 
adequately explained. Further, it is necessary to understand the 
dynamic coupling between the flow and the development of the bottom 
boundary layer before accurate modeling techniques can be developed. 


It is towards this understanding that this work was directed. 


B. PREVIOUS RELATED RESEARCH 

Work by Welander (1957) dealt with obtaining time dependent 
solutions using Ekman dynamics. His results were presented in 
terms of a time dependent Green's function and were related to the 
problems of predicting tides and storm surges. He demonstrated 
that the velocity profile, volume transport, sea level elevation, 


and bottom stress could be exactly represented by his 





"integro-differential" solution. It was suggested that this 
solution could be used to improve storm surge and tide predictions. 
A step-by-step numerical scheme was proposed, but the details of 
the scheme were not presented in the literature. 

Galt (1970) covered much of the same material that Welander 
presented but used a different method of solution. The solution 
allowed a more explicit formulation of results which provided more 
insight into the possible transients present and the relevant time 
Scales of each. He discussed the need for information relating to 
time dependent flow in shallow water regions, and suggested that 
the numerical investigation of the interaction between a seiche 
and its associated transient bottom spiral might be useful in de- 
termining how the various transients present are effectively 


coupled. 


C. THESIS OBJECTIVES 

Since the major area of interest was to be the bottom boundary 
layer, the first objective was to formally solve the governing 
equations neglecting wind stress. The solution would therefore not 
allow a surface Ekman spiral to develop and effect the interior 
fiow. It will be showed that the solution can be represented by 
two equations, One integral and one differential. These equations 
are similar in form to the ones developed by Welander and Galt. 
Although both of these authors proposed that these equations could 
be solved numerically using step-by-step procedures, it was not 
initially apparent whether a stable and accurate numerical scheme 
was feasible. Therefore the second objective was to formulate a 


numerical model capable of simultaneously solving these equations. 
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Finally, as suggested by Galt, this model would be applied to 
represent a seiche to examine the relationships which developed 


between the various transients. 


re 





1k 


A. 


THE FORMAL SOLUTION 


GOVERNING EQUATIONS AND ASSUMPTIONS 


The momentum equations considered are essentially those used 


by Ekman. 


small and neglected. 


coefficient of vertical friction are all assumed constant. 





2 
dSu ou 
oe = 
ot yo 
Zz 
oe + fu - k ~ 
2 


Where u and v are velocities in the 


The non-linear and horizontal friction terms are assumed 


The Coriolis parameter, the density, and the 


Thus: 


(1) 


(2) 


x and y-directions, f the 


Coriolis parameter, g the acceleration of gravity, and k the 


constant coefficient of vertical friction. 


€(x,y,t) is the ele- 


vation of the free surface and the z-axis is considered positive 


up. 


The third equation necessary to 


solve for u, v, and § is the 


vertically integrated continuity equation. 


dg , dU, dv | 
ot (Ox oy : 
where 
O 
U = { WieeGiz =; 
-h 


(3) 


(4) 


is the volume transport and h = h(x,y) is the depth of the water. 


§ is assumed small compared to h. 


2 





Neglecting the effects of wind stress on the motion results in 


the following linearized boundary conditions on u and v: 


(ea (5) 
Ce rn ee (6) 


For this development it was assumed that initially a state 


e 
of rest exists: 


(Up-9 =O 3) (gag FO - (7) 


As in Welander's work the following complex notation was 





introduced: 
w =u + iv 
O§ d& 0€ “a 
dn dx” * dy 
Where i = ¥-1. Using this notation eqs. (1), (2), (6), and (7) 
can be expressed as follows: 
Z 
0 0 : 0 
BET Kt ity = = 9 Be (9) 
Z 
Ow 3 
(8) 2g = 0 (10) 
EG = 0 (11) 
(Wee =O. (12) 


‘Equations (9) through (12) define the initial/boundary valued 
problem that must be solved to specify u and v. Although (32) in 
eq- (9) 1s not externally specified it can theoretically be deter- 
mined by integrating eq. (3). This makes it possible to treat it 


as a time dependent forcing term (i.e. 2 (t)). 
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B. METHOD OF SOLUTION 


To solve the problem begin by multiplying eqs. (9) through 


. : 
(12) by aoe and introducing the relation w = wer? t, Then eqs. 


(9) through (12) can be written as: 





( * 2 * : 
SK = = - gett 2S (t) (13) 
> % 
= )s26 > (14) 
¥ 
(W) p29 = 0 (15) 
(w) op =o, (16) 


A solution to eqs. (13) through (16) can be obtained by using 
a time dependent Green's function approach as described below. 

Assume 

% 
wo = 8 (t) Z(2) (17) 

where Zz 4(2) satisfies the boundary conditions and g(t) satisfies 
the initial conditions. Z_(2) is the eigen-function determined by 
applying standard separation of variables to the homogeneous part 


of eq. (13) using eqs. (14) and (16). It is found to be: 


Z (2) = cos(p2) (18) 


where 


p = (2m-1) 


xf 


(19) 
g(t) is determined by substituting eq. (17) into eq. (13) 
requiring that g(t) satisfy eq. (15), and using the orthogonality 


of the cosine function to determine the appropriate constants. 


Thus 
m pt 2 : : ; 
g(t) = 20d \. o kp (t-t') -ift BB (er yat! (20) 


14 





Therefore a solution to eqs. (13) through (16) is 


m it 2 ; : 
(22) cos (p2) {ke (#8) 488 88 (a ryaes (21) 
p 0 on 


c& 


* 2 
w (x,y,Z,t) = =e m=1 


For the general time dependent solution of eq. (9), dividing 
eq. (21) by eee yields: 


ao 


Z 
mas 52,0) = = » 


; (-1)"cos(pz) er Pot!) at? (22) 
p O n 


C. DISCUSSION OF THE FORMAL SOLUTION 
It is apparent at a glance that eq. (22) satisfies the boundary 
and initial conditions. Substituting eq. (22) into eq. (9) yields 


as a residue: 


2 m 
-1 
ey in cos(pz) =-1. (23) 


It is easily verified that the left hand side of eq. (23) is 
simply the cosine expansion of the right hand side, and therefore 
eq. (22), in the limit, 1S an exact solution to the initial 
problem. 

Equation (22) represents the horizontal currents at all depths 
resulting from the pressure gradient produced by the slope of the 
sea surface at a particular point and instant in time. The integrail 
term carries a running summation of the currents resulting from the 
slope in the past history of the flow. These currents are modi- 


fied in time by a rotation and decay represented by the weighting 


function > 
o7 (kp +if)(t-t'). 


This weighting function introduces two interesting time scales 


into the problem. The first is the inertial rotation period, which 


re 





varies as a function of latitude. For mid-latitude locations it is 
on the order of 20 hours. The second time scale is introduced by 


the exponential decay (T 4) represented by 


2 

1 4h 
tT. 2 = —t (24) 
ie kent) 


This is a function of depth and friction, but for typical values 
(k = .0O1 m“/sec, h = 50m) it is on the order of 28 hours for the 
first term in the series to decay to l/fe of its initial value. 
This is the slowest term to decay. 

The motion represented by (22) depends on the surface slope 
and thus new time scales are introduced into the solution. These 
depend on the horizontal dimensions represented and enter through 
the integration of the continuity equation (3). Time scales re- 
lated to set up of the sea surface and seiche periods will govern 
the time rate of change of the surface slope. This change 
initiates the transients present in the geostrophic and bottom 
currents, which show the inertial rotation and exponential decay 


previously noted. 
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III. THE NUMERICAL SOLUTION 


A. FORMULATION OF THE MODEL 
1. Establishment of a Representative Cross Section 

In formulating the numerical scheme it was decided to 
Simplify the model by neglecting the derivatives alone the y-axis. 
To accomplish this without the loss of any essential physics, it 
was decided to allow the seiche to oscillate across a narrow, 
shallow channel of infinite length. This prevents sea surface 
slopes from developing along the length of the channel but is 
otherwise non-restrictive. The simplification further reduces the 
problem to one of looking at a cross-sectional distribution of 
variables in the channel. 

It was further decided to restrict the model by considering 
a cross-section of constant rectangular dimensions. The width and 
depth of the basin considered were fixed at 5000 meters and 30 
meters respectively, and Figure 1 is a pictorial depiction of this 
cross section. 

These simplifications allow eqs. (3) and (22) to be re- 


presented in the model as 


d o 
28 + — = O (25) 


t 2 : 
a +tif)(t-t") SB (tt )yat' (26) 


gM 8 


2 
w(x,z,t) = = 
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BASIN CROSS-SECTION 





a COO nN 


MODEL CROSS-SECTION 





X=0 | X =k 


Figure 1. Pictorial Depiction of Model Basin 
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2. Establishment of Governing Non-dimensional Ratios 


To represent the various time scales and frictional co- 
efficients possible with more general time dependent motion the 
model was scaled to the period (T) of a free oscillating seiche in 
the model basin, and the coefficient of vertical friction determined 
by Ekman's "depth of frictional resistance (D)'". Equation (26) was 


scaled as follows: 


2L 
ae (27) 
{gh 
a 2k 
D = T F-°? (28) 
Let, 
% 
Cet 
* 
CPS Tt 
* 
dt == T <dt* 
% 
x = = (eleven grid points across the channel) 
*% 
do _ TVA dE 
Ox L > 
ox 


% % ¥% % 
where t , t' , § , x , are non-dimensional quantities and A, is 


the free surface elevation at x = 0, t = 0. Substituting eqs. (27) 


through (29) into eq. (26) yields: 





Zoe 
7 = eo RECs eT eo! ee al if)(t-t') T 
h m=l p L 
0 
*% 
x28 ate (30) 
ox 


: 2m- TT 
Recalling that p = a and 0 = eq. (30) may be written 


FAL 
h T 
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2 2 - 
: LOg Bee co a m i oe -{ 2m (2) zee (>) + i| (ee 


COS 


ae <1 (2m-1)m Os 
*% 
we att (31) 
ox ; 


It was determined that the two non dimensional ratios, 
(f£/0) and (D/h), appearing in eq. (30) govern the time/spacial 
scale and frictional dependence of the motion. For discussion and 
computation the first was designated (PCF) and the second (PCD). 
| The first of these ratios represents the relationship 
between the earth's inertial frequency and the frequency of the 


basin in the following manner. 
f = 2Qsin G 


where {2 is the frequency of the earth's rotation and @ is the 
latitude. 


Cer. 
m 


where ES is the resonance frequency of the basin. Since no external 
forcing terms were used in the solution, Ce is also the approximate 


frequency of the motion. Thus: 


£ Q sing 
CG F 7: 
m 


The model was scaled for (0 S$ PCF <1). A value of 
(PCF = 0) represents motion of a time scale which occurs too 
quickly to be affected by Coriolis accelerations, while a value 
of (PCF = 1) represents motion whose time scale is equal to that of 


inertial rotation and subsequently greatly influenced. 





The second term, PCD, is the ratio of the depth of 
frictional influence (Sverdrup, Johnson, and Fleming, 1942) to the 
actual water depth. The coefficient of vertical friction was de- 
termined in the model from eq. (28) as 


a 
, = (PDe ey Xs ey 


2m 


The model was scaled for (O = PCD S$ 2). A value of (PCD = 0) 
represents a very deep baSin where frictional effects on the motion 
are minimal, while (PCD = 2) represents a basin where friction domi- 
nates the motion. 

If the time scale of the motion being represented by the 
model is controlled strictly by the resonance of the basin, then 
the spacial scale of the basin is determined by a choice of PCD 
and PCF. PCD will govern the approximate depth of the basin while 
PCF determines the baSin width through eq. (27). Certain con- 
clusions about motion caused by external forces (i.e. tides, 
atmospheric pressure) producing time scales independent of the 
basin size were also inferred from the model and are discussed in 
the conclusions. 

3. Adaptation of the Formal Solution to the Model 

To establish two equations that could be solved simulta- 

neously for sea surface elevation and transport, eq. (26) was 


integrated in the vertical using eqs. (4) and (8). Thus: 


0 18 Z 
es, 1 o i = a 
w(x,t) = “SP Bay “234° ep TMD CEY SE (eat a 


and W = U+ iV. Equation (32) represents the total volume trans- 


port associated with the flow at a particular point in the channel. 


Ze 


The real part of eq. (32) represents the cross-channel transport, 
while the imaginary part the along-channel transport. 

The numerical scheme was established to solve eqs. (25) 
and (32) for § and W. The velocity profile, bottom stress, and 
energy balance associated with the flow were also recoverable from 
the scheme, and the details for this recovery will be discussed 
later. 

4. Initial Conditions 
The model assumes an initial state of rest in the motion 


and an initial set-up of the free surface. This was represented 


as 

W(x,0) = 0, (33) 
and 

§(x,0) = a.cos XL , (34) 


where W = W(x,t) and § = €(x,t). 
5. Boundary Conditions 
Since the governing equations were solved without reference 
to horizontal boundaries, these conditions had to be built into the 
model. The model therefore assumes no flow across the horizontal 


boundaries, or 


w(0,t) = W(L,t) = O. (35) 


These conditions were satisfied by requiring the surface slope to 
remain zero at these boundaries at all times. Therefore eq. (32) 
is forced to meet the conditions of eq. (35). 
6. The Step-by-step Scheme 
Both Welander and Galt suggested that a step-by-step 


numerical scheme could be used to solve eqs. (25) and (32). Asa 


22 


first solution the staggered central difference scheme suggested 


by Galt was used and is summarized below. 


t = - At Given initial conditions on W. 

t = 0 Given initial conditions on §. 

t = At Calculate W from eq. (32) using 22 at t = 0. 

t = 2dAt Calculate § using an explicit finite-difference 


scheme on eq. (25) and U at t = At. Calculate a new 
Surface slope. 


ZAt a. MeO Les 


t 
The term At is the time step of the explicit scheme. 

This scheme was developed neglecting Coriolis (f = 0) in 
eq. (32). The scheme proved to be stable and accurate subject to 


the following stability criterion: 


At 
\ ——e < : 6 
gh 7 < 1 (36) 
The expression Ax is the spacial step across the channel. It was 


felt however that a more precise solution, not subject to such a 
restricted stability criterion, could be formulated using a totally 
implicit finite-difference scheme. The development of this method 
1s now presented. 
7. The Implicit Finite-difference Scheme 
The implicit scheme used is a modification of a method 
described by Richtmyer and Morton (1967). 


Start by approximating the time rate of change of the 


transport as follows: 


ry (Waa ng 7 CW) 
oo (37) 
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oy 2. {32g ails se ao(kp tif) (t+at-t") OS ent 


ot At h 1 a 0 ox 
t Bons ; 
2 \. e” (KP at) Lot ) enor’ || (38) 


t fai 
“ ( and regrouping terms, 


By expanding the first integral as { 
O t 


eq. (38) may be written as 


2 , on (kp +if)(t-t')/Q-(kp tif) at_)) erat] 


p “0 2 


t+At ae : 
0 (BB a F ctwteintenee Brenan | aan 


Eqs. (25) and (32) can then be represented in the following form: 


O& 1 | /dU OU 
——) = —-—- = (2) « (Be) | (40) 
ot rst 2 | ox Ox ee 
DW cl, co, [/ve dE 
(>>) an rea Mae wee, (33) ‘ (35 | | (41) 
ort At ZAt Ox ‘ ox mee 
where 
-29 5, 1 (= (kp tif) (t-t"), -(kp tif) at_,, 08 
Cl = h ay 2 ‘ e (e -1) Fit date (42) 
ce ttat 2 
y . <2 il -(kp +1if)(t+at-t' 
PS in & 5, eile vat 
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Equations (40) and (41) can be represented in Matrix form as 


follows: 
oO : B 0 
ee 20s [() . (9) eat | (43) 
2 
where, 
oe 
Cl1/At 
A = D1/At 
O 
0 O C2/At 
B= 0 O D2/At 
-l O O 
and, 
* * 
Cl = Real {cl } G2 = Real {c2 } 
* x 
Die — Imag 1G) } > D2 = Imag (G2 } 


(Real { } and Imag { } represent the real and imaginary part of 
the argument in braces). 
The following implicit finite-difference equation was used 


to represent eq. (43): 


gh tl_ 9” 
J J B Hat bel n An 
at MPS | yell SG 51 eae ae 
where @. = ae 
JAX 
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and n (1,2,3,..., maximum number of time steps) 


i 


Jj (CM 25 Saleen J) 
are the grid points in the time and x-directions respectively. 


Letting 


R 

it 
» 
3 
Q. 


Di ="CGmee +O: - 00. + AtA, 


eq. (44) becomes 


ntl n+l ntl 
- ad. ae + ag. =? 45 
jtl j a2 Jj (45) 


The left hand side of eq. (45) contains the unknown values of 
(9) to be calculated, whereas the right hand side a) are all 


Known values. 


A recursion relation of the following form was assumed. 


n+l n+l 


Qe. = EE. 0, + F, 6 
J J Jrtl J (46) 
where J j j 
Glee ee 
Shing j j 
BE B21 eer <2 3 


and 
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If eq. (46) holds for all eee then the model's left hand boundary 


conditions require that 


0 0 0 O 
Een OO and eo (47) 
oO il 0 


1 


+ - : 
Substituting the expression for Or derived from eq. (46) into 


eq. (45) yields: 


ee | ght, | ai Jo} (48) 
j POE, | Gti | Te EL 
Equating the right hand side of eq. (48) to eq. (46) gives the 
following recursion relationship: 
Of 2 
Be SO eres AL (49) 
7 re ran 
D.-O@F. 1 
= eae) = 1, = ° 
Pj” | Te a a (50) 


and I is the ( 3X3 ) identity matrix. 


Using eqs. (47), (49), and (50), values of EB and F were 


calculated inductively in order of increasing j(j = 1,2,3,..., J-l). 
a . . 
Since the value of a 1S given for j = J-l by the model's right 
hand boundary condition (U. = Vv. = 0nandsc, ee ), values of 
J J J-1 J 


+ 
o. Z were then calculated inductively from eq. (46) in order of 


decreasing j(j = J-l1, J-2, ... , 1). This completes the calculation 


of @” 


J 


representing Bs Fas D5 and eq. (46) are presented in appendix A. 


+ 
: (j = 1,2,3, ... , J-l). The details of the mathematics for 
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8. Numerical Treatment of cl. and C2" 
¥ * 

The numerical scheme for representing Cl and C2 is 
extremely lengthy and is merely outlined in this section. For a 
more detailed presentation the reader is referred to appendix B. 

The integration of cl. was performed by dividing the total 
integration time into small time intervals (At) and assuming that 
the surface slope and effects of inertial rotation could be re- 
presented adequately by a constant throughout the interval. These 
values were calculated at the mid-point of each At interval and 
used in the scheme as representative values throughout the entire 
Sub-interval. This permitted the integral to be evaluated analyt- 
ically over each interval with a maximum error on the order of 
these approximations. Using these assumptions cl was numerically 


represented as , 


MAXM 1-1 2 
ci*t = “225 & [5 a.-e" kp ] [BEN kp at .-ifat_), 


2 
"= (p* L kp N=1 
~if(I--N)At,, -kp“(I-1-N) at 
X (e ive ) (51) 
¥] : 
mer I = 2, 3, ««-« 5, and Cl L = 0. The value of the integer I ranges 


from 1 to MAXT (maximum problem time) and IAt represents the elapsed 
time from initial conditions. L takes on integer values from 1 to 
MAXL (maximum number of cross-channel grid points) and Lax is the 
distance from the left boundary. The integer M ranges from 1 to 
MAXM, and determines the maximum number of modes in the series that 
are included in the approximation of cl. It was determined from 
eq. (23) that MAXM = 25 would represent the solution with a maximum 


error of three percent, and this value was used throughout all 
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numerical evaluations. The integer N takes on the values from 1 
to (I-1) and determines, through the exponential terms, how much 
rotation and decay has taken place in the recent time history of 
each transient considered. Since the terms of the series decay 
more rapidly with increasing friction (especially the higher modes), 
it was not necessary to carry their integration as far back as 
N= 1. It was decided to include only the terms of significant 
magnitude, and this was accomplished by requiring that when the 
value of kp (I-1-N) At in eq. (51) was greater than 10, the cor- 
responding terms were discarded. 

The integration of C2" involves only one time interval 
(t to t+At). Again the amount of inertial rotation was evaluated 
at the mid-point and assumed constant throughout the interval, and 


* 
C2 was evaluated as the following constant: 


z MAXM = ey tn? 
eae _ 22 » “i 1f(sAt) 5 (aree kp ae (52) 
M=1 kp 


¥ ¥ 
Cl and C2 pboth contain the term 


1 -kp~ At ; 

Feat ) » and for k = O this term can 

kp 

not be numerically evaluated in this form. However this difficulty 
-kp“ At 

was Overcome by expanding e in a Taylor series, i.e. 


(1 - kpfar + Lepraty”_ (ip%at)?) (53) 


and substituting this series for the term in eqs. (51) and (52) 


whenever the value of kp“ At was less than 1074, 
The method for separating the real and imaginary parts of 
* * e 3S 
Cl and C2 involves using the relation 


-~if 
ers * = cos(ft) - i sin(ft) 
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in eqs. (51) and (52) and gathering terms. The details are 
covered in appendix B. 
9. Velocity Calculations 
The velocity profiles across the channel were calculated 
by the model at discrete time intervals determined by an input 
parameter. Noting the simularity between the integral of eq. (26) 
and el. it was numerically possible to compute and store the terms 
of these integrals needed for velocity, during the calculation of 
cl. For velocity profile computation, these stored values were 
recalled and applied in eq. (26). The details for numerically re- 
presenting eq. (26) is also presented in appendix B. 
10. Bottom Stress Calculations 
The bottom stress (Tf) associated with the motion was 


calculated using the following relation: 


(mM 
Ss « (3) z=-h 


or from eq. (26), 


© at Dee . 
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Since the integral term of eq. (54) is identical to that 
of eq. (26), bottom stress is easily calculated whenever the 
velocity profile is recovered. Appendix B also contains the 
details for this calculation. 

li. Energy Calculations 
| The total energy balance associated with the seiche is a 


sum of its potential and kinetic energy. This balance was 
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calculated as follows: 


Gb 
PE sp 9 \ e~ ax (55) 
O 


sp 


PE and KE are the potential and kinetic energies and (pf) is the 


0 


a 
\ w dx dz (56) 
O 


KE 
-h 


is the density of the water. Equations (55) and (56) were evaluated 
using Simpson's 1/3 rule (McCalla, 1967) whenever velocity was 
calculated. 

The energy calculations were monitored during model 
testing as a calculational check since small numerical errors 
were quickly apparent in increasing energy levels. 

12. Scaling of Model Output 

Since the governing equations for the model were scaled 
to produce the ratios PCF and PCD, the output of the model must 
be appropriately re-scaled to obtain real values for the variables. 
The series of graphs were produced on a CALCOMP plotter with cor- 
responding values from the model. The various axes were auto- 
matically scaled by the computer to fit the data being represented. 
These axes were then re-scaled using a non-dimensional scaling 
factor derived from the continuity equation to allow real values 
to be easily obtained. The exception is the figures depicting the 
energy balance, which are not scaled and should only be used to 


indicate relative changes in energy levels. 
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The scales for the various axes are as follows: 





AXIS SCALE 


Surface Elevation S/A, 


1. 
2 
Time t (gh) * 


2L 


% 
: 2w_(h/g) 
Velocity / 


CyA, 


2W 
. CA, (gh) 


We 


Transport 


| 


ot (n)2/? 


CA kg@ 
3°09 


Bottom Stress 


The values of g, h, k, A? L are determined from the 


actual basin being represented and C,, C C., are scaling con- 


2 a 
stants dependent on the parameters of the prototype basin and the 


1 


scaling of the graphs by the computer. These values are calculated 
for each figure and inserted in the scale. The variables, t, w, W, 
§, and Ff are the real values associated with the basin being re- 
presented and may easily be obtained from the graphs by entering 
appropriate arguments for g, k, h, Ao: and L in the corresponding 


scale. 


B. TESTING THE MODEL , 


1. Establishing a Reference Solution 
One of the major difficulties in numerical analysis is 
verifying the accuracy and stability of the scheme. Since many 


numerical solutions deal with problems that have no analytical 
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solutions it is common practice to test the accuracy of these 
models against empirical data, or analytical solutions involving 
degenerate cases. As previously noted, data pertaining to time 
dependent flow is extremely scarce, and it was decided to test the 
accuracy of the model against a known analytical solution for the 
degenerate case where PCD = O (i.e. no friction). Analytical 
solutions to this case are readily available in the literature and 
a modification of one proposed by Prodman (1952) was used in the 
verification of the accuracy of the numerical scheme. 
2. Accuracy and Stability of the Model 

The accuracy of the model was determined by comparing the 
numerical representation of the free surface elevation at (x = 0), 
the mid-channel transport, and the total energy balance, with the 
analytical values. Figures 2 through 5 highlite this comparison 
and indicate that the accuracy of the model was a function of time 
step (At) and the Coriolis parameter (f). It was noted that as the 
time step was reduced the numerical solution appeared to converge 
to the analytical solution. The phase shift in surface elevation 
over five periods was reduced from 45 degrees for (At = T/20) to 
24 degrees for (At = T/40) and appeared to be independent of the 
value of Coriolis. This error might possibly be introduced into 
the system by approximating the sea surface slope as a constant 
value in the time history integration scheme. As (At) increases 
this linear approximation is extended over a greater time span 
and the resulting errors amplified. As the value of the Coriolis 
parameter was increased, the error in representing the amplitude 


of the free surface and total energy balance also increased, as 


Se. 
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depicted in Figures 3 and 5. The amplitude deviation was a maximum 
of one percentage for (At = T/20) and the energy rise was 8.5 and 
3.7 percentage for (At = T/20) and At = T/40) respectively. Errors 
of this nature did not appear for solutions which neglected Coriolis, 
and possibly were introduced by approximating the amount of rotation 
ipa . see ail as a constant in the integration scheme. It is 

not exactly clear how this error increases the energy in the model 
but it was observed that it appeared linear and very systematic. 

To determine whether this energy rise constituted a 
numerical instability, a standard Fourier series method (Smith, 
1965) of stability analysis was attempted on eq. (44). However, 
due to the complicated nature of the amplification matrix, the 
eigenvalues could not be readily obtained and the analysis was 
abandoned. As an alternate test, the model was run with (At = T/5) 
for (t = 5T). This corresponds to [fan S¢ > 3 | , and resulted ina 
‘further linear rise in the energy level. These errors were there- 
fore judged to constitute a problem in numerical accuracy and not 
Stability, for the range in which the test cases were run. 

It was further found that the magnitudes of these errors 
could be predicted and controlled by varying the time step, and the 
model was run maintaining a maximum allowable error of one percent 
in the elevation and transport amplitude, and four percent rise in 
the total energy balance. This was accomplished by using a time 
Step of (At = T/20) for (PCF < .5), and (at = T/40) for (PCF 2 .5). 
Since the errors in phase shift were strictly a function of time 


Step, they were completely predictable and accounted for during 


data analysis. 


38 


It should be noted that for cases where friction is 
neglected (k = 0), the transients do not decay with time and must 
be carried throughout the integration scheme. Additionally, the 
bottom boundary condition degenerates to a free slip situation and 
the approximation of a constant velocity profile with a truncated 
cosine series becomes more inaccurate. The number of calculations 
increase rapidly as friction is decreased, which may introduce 
increasing errors in round-off. 

Based on the comparison with the analytical solution it 
was concluded that the accuracy of the model could be sufficiently 
controlled such that the details of time dependent motion, as 
described by the governing integral and differential equations, 


could be adequately represented and examined. 


C. MODEL LIMITATIONS 

The most serious limitation imposed on the model were the 
basic assumptions of Ekman dynamics used in the formal solution 
of the governing equations. These assumptions were considered 
necessary to simplify the mathematical treatment of these equations 
to the extent that a numerical solution was feasible. The additioanl 
Simplifications imposed during the formulation of the model further 
limit its capabilities. It was felt however that valuable experi- 
ence in the techniques of numerically representing time dependent 
motion would be gained during the formulation of this simplified 
model, and insight gained from the numerical product. With this 
information available it might be possible to further refine the 


model and remove many of the more important restrictions. The 


ay 





limitations inherent in the present model are listed in this 
section and the more important ones will be discussed in section 
(V) with recommendations for improving the scheme. 

The basic assumptions made in arriving at eq. (22) have the 
following limiting effects of the model's capabilities: 1) the 
non-linear terms in the motion have been ignored; 2) a variation 
in the coefficient of vertical friction in the layer was not 
permitted; 3) the assumption of constant density neglects the 
effects of stratification on the pressure gradients in the layer, 
and the subsequent modification of the vertical momentum transfer ; 
and 4) horizontal friction was considered insignificant and 
neglected. 

Simplification of the model resulted in the additional 
limitation: 5) the assumption of a constant basin neglects the 
important effects of bottom topography on the flow. This is es- 
pecially important for considering motion in shallow near shore 


areas where these effects are large. 


LO 


IV. PRESENTATION OF RESULTS 


A. TEST CASES AND GEOPHYSICAL SIMULATIONS 
1. Description of Test Cases 
In order to gain maximum insight into the physical processes 
occurring in the motion, as well as to study the effects of inertial 
rotation and friction on the variables, the model was run for 25 
test cases in the anes of (O = PCF = 1) and (0 = PCD S& 2). 
Table I summarized the various controlling parameters of each run 
and includes the Central Process Unit (CPU) time used on the IBM 
360/67 system for each program. 
2. Description of Geophysical Simulations 
As interesting dynamic interactions developed in the test 
cases, it was decided to simulate actual geophysical embayments 
to determine to what extent these phenomena were present and 
Significant. The actual length, depth, and latitude of the basins 
were measured and used to calculate (f), (09), and (PCF). An Ekman 
depth of 20 meters was assumed in the calculation of (PCD). Runs 
26-30 represent these simulations and the various imput parameters 


are summarized in Table II. 


B. DISCUSSION OF TEST RUN RESULTS 
1. Frictional Effects on the Flow 
To describe the conditions where friction totally dominates 
the motion, the model was run with varying amounts of friction and 


(PCF = 0). It must be noted that by neglecting the effects of 


inertial rotation, the Ekman depth is no longer specified as 
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Ay 





Table I. Parameters of Test Cases (runs 1-25) 
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evidenced by eq. (28). Therefore, the frictional coefficients 
generated by runs 1 through 4 were artificially imposed on the 
model with (PCF = 0), and resulted in runs 5 through 8. 

The results for these cases are depicted in Figures 6 
through 8 and are simply a seiche damped by friction. As the 
friction coefficient increased the damping of the motion increased 
with a resulting increase in the period of oscillation. For the 
maximum case of PCD = 2), the seiche was completely damped after 
43; periods. These results were completely predictable and are 
presented to show that the model's representation of motion which 
is dominated by friction conforms to previous observations. 
(Pierson and Neuman, 1966) 

An interesting coupling between bottom stress and transport 
developed during these cases and a typical example is depicted in 
Figure 9. In all cases examined a phase lag quickly developed be- 
tween the transport and the bottom stress. For the case showed, 
the value of this lag averaged 45 degrees through the periods of 
oscillation. This phase lag can be explained by recalling that 
the transport is an integration of the motion over the entire water 
cOlumn, whereas the bottom stress is solely dependent on the 
velocity gradient at the bottom. Since friction has its greatest 
effect at the bottom, changes in magnitude and direction of the 
flow nearly always initially occur next to this boundary and prop- 
agate upwards in the column. The result is that changes in trans- 
port always lag changes in the bottom velocity gradient and 
consequently the bottom stress. This phase lag has obvious effects 


on the accuracy of commonly used numerical schemes where bottom 


Ll 


HO 


an 
oe 


= Add) eoefinsg eeTy oy UO UOTOTIY Fo szOSTTA *g oanbty 


qdod om= ¢ ame w © Gem © awe ¢ a= 


i 


(Oe oe ae 
ddd 





z. 
ae 


“Id 
Y 


4 


OUTIL PoOZT[eWION 
7 ub ya, ai 


T 





O (§/A,) 


Normalized Elevation at x 


45 


O es a7e: 


Sod) zat0dsuezl Teuueyo-ptw 
qgdoOd quae ¢ eee ¢ GEE ¢ eames ¢ a= 


as a a 





dod = 
b ¢ 
. / \ 
~ \ -t “aS: \ 
Sa or a \ 
pee 4% i \\ 
\Vy * 
\ \ 


UO UOTLOTAY Fo 


c au 
a 
\ 
ea 
ax / \\ 
\ : 
yy \ J \ 
‘ / 
fl 2 i| \ 
‘] \ i 
ote x ‘I 
/ : ‘ ye 
Noe 7 hi 
\ / 
NS 


S}JOSOTIA 


ac 


°/ arn6ty 


ses OUTIL pezITeurz0Nn 
} 


a 


\e 


ete anaes 


J 


eo) 


4) W 


~ 
® 


5 
A,(gh) 


Lé 


Normalized Transport 


Z 


re 


"(9 = AOd) aouetTeg A6x9euq UO UOTOT.AY Jo szyOeTJA °Q aanbty 


ic 


SUTL PoOZITeuxoNn 
tio ee 


Cc T O 
O 
oe tan me 
‘. 
Se 
~ 
= 
~ ¢ 
SA 
‘ 
ae ee \ 
a N 
\ 
[a \. 
‘\ \, 9 
~ N, 
SS \. 
~ N 
~ ‘ 
> \ 
\ \ 
\ \ 
‘ \ 
AS \ 6 
NS ‘ 
\ ‘ 
‘\ \ 
~ ‘ 
Sun 
NOS ZT 


a 


a 


Model Energy Balance (joules X 10 


47 


°(T=d0d ‘Q=d0d) }10dsueirL yTauueyo-ptIwn pue ssai3zS woi}0g Fo Uostizeduod 





“Tab ys SUTL pezITewio0N 
E 
G 7 € eC 
I - 
(1) SS3IYS mee ee 
~ I 
By (n) zz0dsuezy i\ 
\ 
AA8, “A ji : 
- g ae /~\ \ ‘ 
f \ / , 
t \, / \ / : / \ 
W . / \ / \ / \ 
“ e / \ ° \ : \ 
- \ / \ ! \ / \ 
» O \ / \ / , : ’ 
‘ \ / \ / \ / \ 
5 \ Zo a \ / \ i \ 
is : > ~, / \ / \ 
a : VV j ’ 
TO . \ / \ 
Y \e 
2 
2 \ 
oO 
E 
H 
Oo 
es 
i 


°6& ainbty 


2 


5.4) W 
A, (gh) 


Normalized Transport 


- 48 


stress is parameterized in terms of the transport. This technique 
places strong restrictions on possible interactions between the 

flow and the bottom, and although it may be sufficient for approx- 
imating steady state problems, it is unlikely that it will represent 
transient motion very well. 

Figure 10 represents typical velocity profiles for motion 
which is totally dominated by the effects of vertical friction. 

The increased period of oscillation and the mechanism causing the 
phase lag between the transport and the bottom stress is clearly 
evident. At times (t = T/2) and (t = 3T/4) the profiles are in 
transition from being totally positive to totally negative, and the 
transport and bottom stress are 180 degrees out of phase. 

These velocity profiles are typical of those found in other 
boundary layer flow and are similar to the wind profiles observed 
at the air sea interface over smooth water. They indicate that the 
stability of the flow has a strong frictional dependence. At 
(t = T/4) and (t = T) the profiles are maximum and, except initially 
in the upper layers, friction exerts a strong influence throughout 
the column. The result is a smooth regular flow and a well de- 
veloped frictional boundary layer with the region of greatest shear 
at the bottom. During the transition phases (t = T/2 to t = 3T/4) 
the profile changes direction and the water column appears to be 
less stable. The velocities in the column are minimum just prior 
to and after reversal and thus the influence of friction in the 
column is reduced. This region of reduced frictional influence 
moves up the column at the point where the flow changes direction. 


As the flow completes its reversal and the velocities in the column 
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Figure 10. Effects of Friction on Mid-channel Velocity 
(PCF = Of PCD = 2). 
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increase, the effect of strong frictional influences is quickly 
re-established throughout the layer, the region of maximum shear 
shifts to the bottom, and smooth flow once again prevails at time 
(t =T). 

2. Coriolis Effects on the Flow 

To gain insight into how inertial rotation modifies the 
transient motion, runs 10, and 23 through 25 were made neglecting 
friction (PCD = 0). 

It is important to realize that the relationship between 
the transport and the free surface elevation for a non rotational, 
frictionless seiche, as seen in Figures 6 and 7, is similar to the 
velocity/elevation relationship of the frictionless pendulum. At 
(t = O), (t = T/2), and (t = T), the energy of the seiche is en- 
tirely potential, whereas at (t = T/4) and (t = 3T/4) it is 
rely kinetic. At intermediate times the energy balance is a 
combination of potential and kinetic and is totally conserved. The 
free surface elevation and transport periodically oscillate between 
maximum positive and negative values of respectively equal magni- 
tudes, always maintaining their constant phase shift. Figure 2 
represents this oscillation for the free surface. 

As the effects of Coriolis acceleration were introduced 
into the model these conservative oscillations were greatly modified 
as can be seen from a comparison of Figures 2, 3, 11, and 12. The 
effect of a moving earth is that the total potential energy of the 
Surface elevation will not be available for transfer to kinetic 
energy Since the water particles are forced to turn due to the 


Coriolis Acceleration. Although the energy balance still remains 
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Figure ll. Effects of Inertial Rotation on Surface Elevation (PCD = 0). 
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conservative, the relationship between the free surface elevation 
and the flow changes Significantly. 

To understand how rotation acts to modify the flow, con- 
sider first the physical processes which drive the system. These 
are discussed in terms of the frictionless irrotational Seiche. 

The initial conditions for the model were identical for all runs; 
that of a positive set-up of the free surface at the left boundary 
(x - 0) and no motion in the basin (a condition of static equili- 
brium). As the system is released the pressure gradients developed 
by this set-up start the fluid moving in the positive x-direction 
resulting ina net divergence in the left side of the channel and 
convergence in the right. The free surface responds to this flow 
by slumping at the left and rising at the right with equal magni- 
tudes. At time (t = T/4) the surface is level (i.e. no gradients) 
and the kinetic energy associated with the flow is equal to the 
potential energy of the initial free surface set-up. Therefore the 
surface continues to fall at the left and rise at the right until 
this kinetic energy has been entirely converted into potential 
energy and a state of momentary equilibrium is once again reached. 
The sea surface is now set up at the right boundary equal in magni- 
tude to the initial conditions, the pressure gradients equal but 
opposite to the initial values, and the flow reverses direction 

and the system returns to initial conditions, etc. It is important 
to realize that only the divergence of the flow in the cross-channel 
direction can influence the system because no gradients are allowed 


to develop along the length of the channel. 
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With these processes in mind the observed effects of ro- 
tation on the system can adequately be explained. With the initial 
Set-up previously described, the flow will start in the same manner. 
As soon as the motion begins however, the Coriolis acceleration will 
start to turn the flow in the negative y-direction (left rotation), 
and a component of the flow will start to develop along the axis of 
the channel. Therefore some of the potential energy is converted 
into rotational kinetic energy. The amount of rotational kinetic 
energy that the system receives is a function of the velocity of 
the flow, the value of Coriolis, and the time scale of the motion. 
The slumping of the surface elevation at the left and rising at the 
Yight will continue as long as there is divergence in the positive 
x-direction; however since a certain amount of energy is rotational, 
the set-up at the right boundary can not achieve the magnitude of 
initial conditions. When the cross-channel flow vanishes the en- 
tire motion is in the along channel direction, with subsequent 
rotation turning the flow back in the direction from which it 
initially came. It should be noted that for a maximum value of 
(PCF = 1), Figure 3 shows the free surface elevation at the left 
boundary never falling below the horizontal. 

An interesting relationship developed in the transport as 
rotation was introduced into the system, resulting in a rectifi- 
cation of the along channel component of the transport (V). An 
examination of Figure 12 showed that the cross-channel transport (U) 
periodically oscillates between maximum positive and negative values 
while the along channel component oscillated between zero anda 


maximum in one direction. It was easily demonstrated that this 
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direction is strictly a function of initial set-up. As the flow 
was initiated in the positive x-direction the initial along channel 
component developed in the negative y-direction. When the cross- 
channel flow reverses direction, along channel components are 
generated in the positive y-direction. These components are just 
sufficient to cancel the components previously established in the 
opposite direction, resulting in a rectification in the along- 
channel transport. For motions of large time scales, this process 
implies that important changes in the flow are caused by Coriolis 
accelerations and these will be discussed in the etna tne 

Another observed feature of rotation was the subsequent 
reduction in oscillation period for increased values of Coriolis. 
The effect of rotation in reducing the spacial scale of the motion 
allows the system to oscillate with greater frequency. 

It 1S apparent from the previous discussions that the 
effects of inertial rotation can not be neglected when dealing with 


large time scale motion. 


3. Combined Effects of Friction and Coriolis 


The remainder of the test runs were used to investigate 
how friction and inertial rotation in various combinations effect 
the transient motion. Some of the important consequences are de- 
picted in Figures 13 through 17. 

The phase shift difference between bottom stress and 
transport existed throughout all values of (PCF) and (PCD), as 
evidenced in Figures 13 and 14. This phase lag occurred in varying 
degrees but was a definite characteristic of time dependent motion, 


vanishing only when the motion became critically damped. 
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The rectification of the along channel flow was also a pre- 
dominate characteristic of the flow where inertial rotation was 
considered. Figures 13 and 14 again show this relationship, as 
modified by friction, and Figure 15 depicts the combined effect on 
the free surface. It can be seen that for an extreme case (PCD = l, 
PCF = 1) the free surface only very slowly slumps to a zero potential 
because friction retards the cross-channel flow long enough for ro- 
tation to reverse its direction before the free surface reaches the 
horizontal. 

Figure ieee that with rotation the energy level is 
damped with relatively fewer oscillations. It must be realized 
however that as (PCF) increases, the model is describing motion on 
an increasing time scale. 

Figure 17 represents the path that would be taken by fluid 
particles (or suspended matter) at three levels in the column over 
five periods of oscillation. In addition to the rectification pre- 
viously noted, there appears to be a net drift of the motion in the 
initial cross-channel transport direction. A closer examination of 
Figures 13 and 14 also show this drift and is a combined result of 
rotation and friction. Since the motion is being continuously 
damped by friction, the cross-channel transport can not return the 
fluid to its initial position after each oscillation and a net 
transport in the initial direction of motion results. The net 
drift is then strictly a function of initial set-up conditions of 
the free surface. The velocity profiles depicted in Figures 18 
through 20 show some of the fine details of the flow with the 


combined influences of friction and inertial rotation. 
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Figure 18. Effects of Friction and Inertial 
Rotation on Mid-channel Velocity 
(BCR = sree eCbe- at). 
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Figure 19. Effects of Friction and Inertial Rotation 
.. On Mid-channel Velocity (PCF = 1, PCD = 1, 
t = 3T/4). 
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Figure 20. Effects of Friction and Inertial Rotation on 
Mid-channel Velocity (PCF = 1, PCD =1, t = T). 
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The controlling influence on the stability of the water 
column continues to be the friction, however the interactions of 
the fluid are now more complicated due to the introduction of 
inertial rotation. Under the influence of increasing frictional 
dominance towards the bottom, the left rotation induced in the flow 
by Coriolis is increasingly damped, resulting in the familiar Ekman 
spiral. During periods of the motion previously described as 
smooth, the spiral remains fairly constant as the column rotates 
to the left. However, during the transition periods (Figure 19) 
when frictional influence is reduced due to decreasing flow velocity, 
the spiral unwinds from the bottom and varing amounts of twisting of 
the column results. This is again evidence of the transitional 
nature of the fluid during these periods. It should be noted that 
the velocity scale of Figure 19 is an order of magnitude smaller 
than that of Figure 18, indicating how small the velocities of the 
colLumn are during transition. As the flow once again becomes uni- 
form, the spiral is quickly re-established under the stabilizing 
influence of bottom friction, as depicted in Figure 20. 

It is interesting to note from Figure 18 that during the 
transition periods the along channel component of velocity does in 
fact change direction for brief periods. Because these periods are 
of such short duration, this velocity reversal along the channel does 


not contribute to any significant transport, in this direction. 


C. DISCUSSION OF THE GEOPHYSICAL SIMULATIONS 
All of the previously discussed processes characteristic of 


time dependent motion, were observed during runs 26 through 30. 
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However, Since the magnitude of these processes are Coriolis and 
frictional dependent, they occurred on a much smaller scale than 
observed for the test cases. 

The model generated the time dependent variables produced by 
the seiche on a time scale equivalent to the resonance period of 
the eit For these geophysical simulations, the periods were 
all too small to allow inertial rotation to have any significant 
effect on the motion. However, for shallow seas of large hori- 
zontal extent, values of (PCF) approaching .25 are entirely pos- 
sible and test runs 1 through 4 showed that inertial effects would 
be significant and must be considered. Additionally, motion whose 
time scale is driven by external forces, independent of the re- 
sonance basin, are likely to be significantly altered by rotational 
dependence. For example the long period waves propagating up and 
down the continental shelf probably experience all the rotational 
characteristics associated with the time dependent motion described 
in the previous section. The observed phenomenon of transport 
rectification might occur in these systems and suggests a mechanism 
for modifying bed load transfer across the shelf. This process is 
a function of initial conditions and appears to require a period of 
static equilibrium just prior to initial motion. A region where the 
seasonal migration of atmospheric pressure systems are predominantly 
in One direction, might satisfy these requirements, with a resulting 
net rectification in One direction. Although this process is highly 
speculative, it certainly appears possible and a modification to 
the model to facilitate its study is suggested in Section V. 


The frictional effects observed during runs 26 through 30 were 


also minimal. It should be noted that Ekman depths of greater than 
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20 meters are possible and thus these test runs represent a minimum 
frictional dependence. 

The exception to these minimal effects was the observed phase 
lag between transport and bottom stress. This lag was pronounced 
in all cases and for the simulation of South San Francisco Bay, 
for example, averaged 40 degrees. This indicates that even small 
amounts of friction are important and re~-emphasizes the need for a 
mathematical representation of time dependent transients that allows 
independent calculation of variables such as bottom stress and 


transport. 
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V. CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER WORK 


The model demonstrated that eq. (22) was a suitable mathematical 
representation of the time dependent motion described by the governing 
equations and the initial/boundary conditions. In this formulation 
the time dependent variables were uniquely determined by the local 
time histories of the surface at each point and this allowed inde- 
pendent representation of the various transients. The model further 
demonstrated the feasibility of developing a totally implicit finite- 
difference scheme for approximating eqs. (3) and (22) with accept- 
able accuracy. With the insight gained from the results it should 
be possible to improve the accuracy of the system by removing some 
of the important restrictions imposed by the Ekman dynamics and 
the simplified model. 

The physical processes observed with time dependent motion 
showed the relative importance of friction and inertial rotation 
in determining the nature of the variables. It was apparent that 
for the larger time scales (PCF on the order of .25 or larger), 
Coriolis accelerations are important to the flow and must be con- 
Sldered. However, shorter time scales will probably be adequately 
represented neglecting rotation with little accuracy loss. During 
the data analysis it appeared that friction, even in small amounts, 
exerted considerable influence on the characteristics of the vari- 
ables and their relationships, and should therefore not be neglected 


in representing time dependent motion. 
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Since frictional effects are predominant in the system, the 
simplified assumption of a constant coefficient of vertical friction 
throughout the layer, is probably one of the most limiting factors 
on the model's accuracy. A technique for improving the repre- 
sentation of the influence of friction might be the assumption of 
constant stress layer (Prandl layer) underlying the Ekman layer. 
The solutions in these layers would be matched at their interface; 
that 1S, at some height above the bottom determined as a function 
of the amount of friction present and the roughness of the bottom 
boundary. This method has been used in air-sea boundary layer 
models with considerable improvement in accuracy, and merits con- 
Sideration at the ocean boundary in this scheme. 

The assumption of constant depth does not account for modifi- 
cations of the flow caused by changes in the potential vorticity 
induced by the shrinking and stretching of the water column over 
varying topography. This restriction is probably most serious in 
shallow regions where time dependent flow is significant and changes 
in the bottom extreme, for example the flow in the vicinity of the 
Monterey Submarine Canyon. 

Certain large scale motions merit representation and investi- 
gation using this formulation. One of these is the study of 
seiches, especially the trapped modes, on the continental shelf 
Since their buildup can cause storm surge (Pierson and Neuman, 
1966). Also, as previously mentioned, the large scale oscillations 
on the continental shelf might make important contributions to 
sediment transfer across the shelf. Therefore it might prove use- 


ful to remove the horizontal boundary conditions from the model 


» 
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and drive the flow by an externally specified time dependent 
forcing term. This would indicate how effectively the motion was 
coupled to these terms and show whether the previously noted trans- 
port rectification is of sufficient magnitude to be significant. 

An examination of Tables I and II indicates how costly this 
type of numerical representation is in terms of computer time. It 
was not determined how much loss in resolution would occur by in- 
creasing the spacial grid size and reducing the length of the local 
time history considered. It is clear that considerable reduction 
in computing time would result and that such considerations are 
worth investigating before a similar model is used for extensive 


production runs. 
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APPENDIX A 


DETAILS OF THE FINITE-DIFFERENCE SCHEME 


Consider first the calculation of the components of the vector 


Di? recalling that 
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Expanding eq. (Al) using (A2) through (A4) yields: 
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Recall eq. (49) and express it in inverse form as 
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where I is the (3 X 3) identity matrix. Using the appropriate 


values for I, @, and Sea it is easily showed that 
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(I + OB; 1) = Zi aps(y27+) (All) 


Where pet (v2 ~*) is the determinant of the matrix (yJ7t) and 
ADJ (y2 7+) is its adjoint. The computation of eq. (All) are ex- 
tremely lengthy and are not included. The resulting product of 


applying eq. (All) to eq. (A9) is as follows: 
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This completes the calculations of the components of Ey: 


Finally consider the calculations of the components of the vector 
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F.= § £3 (A14) 


Recall eq. (50) and express in its inverse form as 
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Again the computations of eq. (All) are involved and are not 


included. The results of applying eq. (All) to (A15) are as 
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Note that the middle term of the vector Be is the entire quantity 
in the brackets. 
The implicit scheme for computing U, V, and § can now be 


constructed. Writing eq. (46) as 


nt+lL n+l 
G. = EE. _O. + F, Al8 
aid g-L gy j-1 ( 


and substituting the quantities previously derived for 5 1 and 


Pj-1 yields: 
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Equations (A7), (Al2), (A113), (Al6), (A17), and (A19) through 
(A21) represent the expressions used to implicitly calculate values 


n , . : ; : 
for ae at each grid point in the manner described in the text. 
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APPENDIX B 
se #% 
DETAILS OF REPRESENTING Cl _ , C2 , VELOCITY AND BOTTOM STRESS 


¥% 
Consider the representation of Cl by first noting that the 


limits of integration have the following meaning: 


t = O represents the initial condition time. 
t = t represents the present time. 
t = t + At represents the time at which the value of the 


various variables are to be calculated. 

The integration of eq. (48) was performed numerically by 
dividing the total time into small intervals (At) and evaluating 
P(t") as a constant at the mid-point of each interval. This 


5 
allowed Cl to be numerically represented as follows: 
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With the further assumption that the effect of inertial ro- 


: : : -if(t-t' : 
tation was constant in each interval, e ( was approximated 


at the mid-point as eee cat and treated as a constant. 
Analytical integration of the remaining integral resulted in 

eq. (51). To solve for Cl and Dl the following relation was used, 
and the real and imaginary terms collected. 


elft = cos(ft) - i sin(ft) (B2) 
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This allowed Cl and Dl to be represented as 
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Using the same assumption regarding constant inertial rotation 
¥% 
in each time interval, C2 readily reduced to eq. (52). Applying 


eq. (B2) to eq. (52) yields: 
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Consider now the scheme formulated to recover the velocity 
profile. Since velocity was calculated at time t + At, eq. (26) 


was represented by: 


MAXM ,_,,m a eee _ guys ee mnie 
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Applying the same assumptions and relations to eq. (B/7) as 


were applied to eq. (Bl) yields: 
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The term (A5) represents the depth dependence and (A6) the 
slope in the most recent time interval. (A6) was calculated and 
stored by the model after U, V, and § were computed. 

The final quantity needing representation is bottom stress (T). 


Recalling eq. (54) and noting that the integral is identical to 


Ve 


that of eq. (26), bottom stress is easily represented as: 
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Where re) and re are cross-channel and along-channel components 
of stress respectively. 

During the calculation of Cl, Dl, C2, D2, the model stores 
those common quantities needed for velocity and stress calculations 
as defined by eqs. (B3) through (B6), and (B8) through (Bll). The 
program recalls these quantities and applies them to eqs. (B8) 
through (Bll) whenever it is required to calculate velocity and 


bottom stress. 
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APPENDIX C 
COMPUTER PROGRAM DESCRIPTION 

The computer program was written in FORTRAN IV and was used 
with the IBM 360/67 computer system at the W.R. Church Computer 
Center, Naval Postgraduate School. The program has been generalized 
such that basin size, run time, initial free surface elevation, and 
number of modes considered in the solution can be varied in the data 
statements. Additionally, the model allows variable grid spacing, 
requiring only changes in the imput parameters and DIMENSION state- 
ments. The influence of inertial rotation and vertical friction are 
introduced into the model by the two non dimensional ratios dis- 
cussed in the text. 

FORTRAN IV symbols used for the major model parameters are 
defined below. 
A Model time (real time normalized to the free oscillating 


seiche period) 


DH Spacial step in depth (fixed at one meter) 
DX Spacial step in across-channel direction 
DI Time step 

ETA Free surface elevation 


ETAX Sea surface slope 


ETAXS Array for storing the time history of the sea surface 


slopes 
F Coriolis parameter 
FREQ Frequency of the free oscillating seiche 
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HI 


MAXH 


MAXL 


PCD 


PCF 


PCT 


ie 


PERD 


RHO 


SE 


SIGMA 


STRESC 


STRESR 


T 


IE 


TOT 


UT 


VELC 


VELR 


VT 


Acceleration of gravity 

Depth of water 

Free surface elevation at left boundary 

Integer representation of water depth (must be even) 
Maximum number of cross-channel grid points (must be odd) 
Interval in time step when velocity, energy, and bottom 
stress is calculated 

Maximum number of time steps 

Ratio of Ekman depth to depth of the water 

Ratio of Coriolis parameter to depth of the water 

Number of time steps per free wave period 

Potential energy of the seiche 

Period of a free oscillating seiche 

Coefficient of vertical friction 

Density of the water 

Kinetic energy associated with the y-component of velocity 
Angular frequency of a free oscillating seiche 

y-component of bottom stress 

x-component of bottom stress 

Real time 

Kinetic energy associated with the x-component of velocity 
Total energy balance of the seiche 

x-component of volume transport 
y-component of velocity 
x-component of 


velocity 


y-component of volume transport 
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The procedure for initialization of the program is as follows: 

1) Modify the initial DATA statements to correspond to the 
parameters of the motion being represented. 

2) Change DIMENSION statements accordingly. 

3) Insure that appropriate storage and computer time require- 
ments are satisfied. The program required a maximum of 100K BYTES 
storage on the IBM 360/67 system and variable CPU usage time as 
depicted in Tables I and II. 

4) A minor modification to the program is necessary when 
frictionless motions are being represented. If (PCD = O) or 
(PCF = 0), remove the following three expressions from the program. 

K = (I-l1) - 10.0/cc 

Thitereo — 1 

po 40 J = kK, IMl 
Replace these expressions by 

po 40 J = 1, IMl 

5) Because of the formulation, the program can not compute 
velocity at I = 1; therefore a value of MAXP = 1 should not be 
used. To recover the velocity profile at each time step (I # 1), 
initialize the DATA with MAXP = 2 and insert the statement 
(MAXP = 1) immediately before statement number 300. 

The program, aS appearing in the next section, is iniriaieens 
in the MKS system and changes in the imput should be made ac- 
cordingly. The major computational sections of the program are 
prefaced by appropriate comment statements, and a reader with 
FORTRAN IV experience should have little difficulty following the 


program. 
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COMPUTER PROGRAM 
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